/*
   Copyright (c) 2011, The Chinese University of Hong Kong

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/

/*
   This software contains source code provided by NVIDIA Corporation.
*/

/**
 * @file CUDAagrep.cpp
 *
 * @brief This is a CUDA implementation of the agrep algorithm,
 * which is used for approximate string matching.
 *
 * @author Hongjian Li, The Chinese University of Hong Kong.
 *
 * Copyright (C) 2011 Hongjian Li, The Chinese University of Hong Kong.
 */

#include <cutil_inline.h>	/** Used for timing and error handling by auxiliary macros such as cutilSafeCall. */

// About corpus.
#define MAX_CORPUS_FILE_COUNT 40	/**< Maximum number of corpus files, i.e. maximum number of lines in Corpus.txt. */
#define MAX_CORPUS_FILE_PATH_LENGTH 100	/**< Maximum length of corpus file path, i.e. maximum length of lines in Corpus.txt. */
char corpus_file[MAX_CORPUS_FILE_COUNT][MAX_CORPUS_FILE_PATH_LENGTH + 1];	/**< Corpus files. */
unsigned int corpus_file_count;		/**< Actual number of corpus files. */

// About query. One query consists of a pattern followed by an edit distance.
#define MAX_PATTERN_COUNT 10000	/**< Maximum number of patterns, i.e. up to this many lines of matches will be stored into Result.csv. */
#define MAX_PATTERN_LENGTH 64	/**< Maximum pattern length. A 64-bit machine word has 64 bits. */
char pattern[MAX_PATTERN_COUNT][MAX_PATTERN_LENGTH + 2];	/**< Patterns. */
unsigned int edit_distance[MAX_PATTERN_COUNT];	/**< Number of primitive operations necessary to convert one string into an exact copy of the other. Primitive operations include insertion, deletion, and substitution of one character. */
unsigned int query_count;	/**< Actual number of queries. */

// About sequence.
#define MAX_SEQUENCE_COUNT 10000/**< Maximum number of sequences in Corpus.txt. Note that one .fa file might contain more than one sequence if it is not an assembled chromosome. */
#define MAX_LINE_LENGTH 1000	/**< Maximum length of lines of fasta files. Fasta files are scanned line by line in this program. This value should be larger than or at least equal to the length of fasta header, which exceeds far 70 in most cases. */
char sequence_header[MAX_SEQUENCE_COUNT][MAX_LINE_LENGTH];	/**< Headers of sequences. */
unsigned int sequence_length[MAX_SEQUENCE_COUNT];	/**< Lengthes of sequences. */
unsigned int sequence_cumulative_length[MAX_SEQUENCE_COUNT + 1];	/**< Cumulative lengths of sequences, i.e. 1) sequence_cumulative_length[0] = 0; 2) sequence_cumulative_length[sequence_index + 1] = sequence_cumulative_length[sequence_index] + sequence_length[sequence_index]; */
unsigned int sequence_count;	/**< Actual number of sequences. */

// About nucleotide.
#define CHARACTER_CARDINALITY 4	/**< One nucleotide is either A, C, G, or T. */
#define MAX_SCODON_COUNT 220000000	/**< Maximum number of special codons. The current value is 0.22G, therefore the special codon array can store up to 3.52G (=0.22G*16) base pairs. */
unsigned int scodon[MAX_SCODON_COUNT];	/**< The entire genomic nucleotides are stored into this array, one element of which, i.e. one 32-bit unsigned int, can store up to 16 nucleotides because one nucleotide can be uniquely represented by two bits since it must be either A, C, G, or T. One unsigned int is called a special codon, or scodon for short, because it is similar to codon, in which three consecutive characters of mRNA determine one amino acid of resulting protein. */
unsigned int character_count;	/**< Actual number of characters. */ 
unsigned int scodon_count;	/**< Actual number of special codons. scodon_count = (character_count + 16 - 1) >> 4; */

// About result.
#define MAX_MATCH_COUNT 10000	/**< Maximum number of matches of one single query. */

// About CUDA implementation.
#define MAX_UNSIGNED_INT	0xffffffffUL	/**< The maximum value of an unsigned int. */
#define MAX_UNSIGNED_LONG_LONG	0xffffffffffffffffULL	/**< The maximum value of an unsigned long long. */
#define B 7	/**< Each thread block consists of 2^B (=1<<B) threads. */
#define L 8	/**< Each thread processes 2^L (=1<<L) special codons plus those in the overlapping zone of two consecutive threads. */
// Since each thread block processes 1 << (L + B) special codons, the number of thread blocks will be up to (MAX_SCODON_COUNT + 1 << (L + B) - 1) >> (L + B).
// This program uses 1D CUDA thread organization, so at most 65,536 threads can be specified.
// Therefore, the inequation ((MAX_SCODON_COUNT + (1 << (L + B)) - 1) >> (L + B)) <= 65,536 must hold.
// MAX_SCODON_COUNT = 0.22G ==> L + B >= 12 is required.
#define MAX_BLOCK_COUNT MAX_SCODON_COUNT >> (L + B)	/**< Maximum number of thread blocks. */
unsigned int block_to_sequence[MAX_BLOCK_COUNT];	/**< Mapping of thread blocks to sequences. */
unsigned int block_count;	/**< Actual number of thread blocks. */

/**
 * Transfer necessary parameters to CUDA constant memory.
 * This agrep kernel initialization should be called only once for searching the same corpus.
 * @param[in] scodon_arg The special codon array.
 * @param[in] character_count_arg Actual number of characters.
 * @param[in] match_arg The match array.
 * @param[in] max_match_count_arg Maximum number of matches of one single query.
 */
extern "C" void initAgrepKernel(const unsigned int *scodon_arg, const unsigned int character_count_arg, const unsigned int *match_arg, const unsigned int max_match_count_arg);

/**
 * Transfer 32-bit mask array and test bit from host to CUDA constant memory.
 * @param[in] mask_array_arg The mask array of a pattern.
 * @param[in] test_bit_arg The test bit.
 */
extern "C" void transferMaskArray32(const unsigned int *mask_array_arg, const unsigned int test_bit_arg);

/**
 * Transfer 64-bit mask array and test bit from host to CUDA constant memory.
 * @param[in] mask_array_arg The mask array of a pattern.
 * @param[in] test_bit_arg The test bit.
 */
extern "C" void transferMaskArray64(const unsigned long long *mask_array_arg, const unsigned long long test_bit_arg);

/**
 * Invoke the cuda implementation of agrep kernel.
 * @param[in] m Pattern length.
 * @param[in] k Edit distance.
 * @param[in] block_count Number of thread blocks.
 */
extern "C" void invokeAgrepKernel(const unsigned int m, const unsigned int k, const unsigned int block_count);

/**
 * Get the number of matches from CUDA constant memory.
 * @param[out] match_count_arg Number of matches.
 */
extern "C" void getMatchCount(unsigned int *match_count_arg);

/**
 * Encode a character to its 2-bit binary representation.
 * The last two but one bits are different for A, C, G, and T respectively.
 * Note that some genomes contain 'N', which will be treated as 'G' in this encoding function.
 *
 * 'A' = 65 = 01000<b>00</b>1
 *
 * 'C' = 67 = 01000<b>01</b>1
 *
 * 'G' = 71 = 01000<b>11</b>1
 *
 * 'N' = 78 = 01001<b>11</b>0
 *
 * 'T' = 84 = 01010<b>10</b>0
 * @param[in] character The character to be encoded.
 * @return The 2-bit binary representation of given character.
 */
inline unsigned int encode(char character)
{
	return (character >> 1) & 3;
}

/**
 * Parse the corpus container file. It is the first argument to the program. By default it is Corpus.txt.
 * @param[in] corpus_container_file The container file specifying all the corpus files.
 */
void parseCorpusContainerFile(const char *corpus_container_file)
{
	FILE *corpus_container_file_handle = fopen(corpus_container_file, "r");	// The handle to corpus container file.
	corpus_file_count = 0;
	while (fgets(corpus_file[corpus_file_count], MAX_CORPUS_FILE_PATH_LENGTH, corpus_container_file_handle) != NULL)
	{
		corpus_file[corpus_file_count][strlen(corpus_file[corpus_file_count])-1] = '\0';
		corpus_file_count++;
	}
	fclose(corpus_container_file_handle);
}

/**
 * Build the special codon array. Corpus files are enumerated. Statistics of the corpus are calculated.
 */
void buildSpecialCodonArray()
{
	FILE *corpus_file_handle;	// The handle to corpus file.
	unsigned int corpus_file_index; // Used to enumerate corpus_file[].
	char line_buffer[MAX_LINE_LENGTH];	// fgets(line_buffer, MAX_LINE_LENGTH, corpus_file_handle);
	unsigned int line_buffer_length;	// Length of line_buffer.
	unsigned int line_buffer_index;	// Used to enumerate line_buffer[].
	int sequence_index;	// Used to enumerate sequence_header[] and sequence_length[].
	unsigned int character_index;	// Index of the current character across all the sequences of the entire corpus.
	unsigned int character_index_lowest_four_bits;	// The lowest four bits of character_index;
	unsigned int scodon_buffer;	// 16 consecutive characters will be accommodated into one 32-bit unsigned int.
	unsigned int scodon_index;	// scodon[scodon_index] = scodon_buffer; In CUDA implementation, special codons need to be properly shuffled in order to satisfy coalesced global memory access.

	sequence_cumulative_length[0] = 0;
	sequence_index = -1;
	character_index = 0;
	scodon_buffer = 0;
	for (corpus_file_index = 0; corpus_file_index < corpus_file_count; corpus_file_index++)
	{
		corpus_file_handle = fopen(corpus_file[corpus_file_index], "r");
		while (fgets(line_buffer, MAX_LINE_LENGTH, corpus_file_handle) != NULL)
		{
			line_buffer_length = strlen(line_buffer) - 1; // Exclude the last character which is '\n' of ASCII code '\10'.
			if (line_buffer[0] == '>') // Header line
			{
				if (++sequence_index > 0) // Not the first sequence.
				{
					sequence_cumulative_length[sequence_index] = character_index;
					sequence_length[sequence_index - 1] = character_index - sequence_cumulative_length[sequence_index - 1];
				}
				memcpy(sequence_header[sequence_index], line_buffer, line_buffer_length);
				sequence_header[sequence_index][line_buffer_length] = '\0';
			}
			else
			{
				// Use line_buffer_index to enumerate line_buffer, equivalent to foreach (char character in line_buffer)
				for (line_buffer_index = 0; line_buffer_index < line_buffer_length; line_buffer_index++)
				{
					character_index_lowest_four_bits = character_index & 15;
					scodon_buffer |= encode(line_buffer[line_buffer_index]) << (character_index_lowest_four_bits << 1); // Earlier characters reside in lower bits, while later characters reside in higher bits.
					if (character_index_lowest_four_bits == 15) // The buffer is full. Flush it to the special codon array. Note that there is no need to clear scodon_buffer.
					{
						scodon_index = character_index >> 4;
						// scodon_index can be splitted into 3 parts:
						// scodon_index = block_index << (L + B) | thread_index << L | scodon_index;
						// because 1) each thread block processes 1 << (L + B) special codons,
						//     and 2) each thread processes 1 << L special codons.
						// The scodon_index is accommodated in the lowest L bits.
						// The thread_index is accommodated in the middle B bits.
						// The block_index  is accommodated in the highest 32 - (L + B) bits.
						// This program uses 1D CUDA thread organization, so at most 65,536 threads can be specified, i.e. at least 16 bits should be reserved for block_index.
						// Therefore, the inequation 32 - (L + B) >= 16 must hold. ==> L + B <= 16.
						// In order to satisfy coalesced global memory access, thread_index should be rearranged to the lowest B bits.
						// To achieve this goal,
						//         1) block_index remains at the highest 32 - (L + B) bits,
						//   while 2) thread_index should be rearranged to the lowest B bits,
						//     and 3) scodon_index should be rearranged to the middle L bits.
						// Finally, scodon_index = block_index << (L + B) | scodon_index << B | thread_index;
						scodon_index = (scodon_index & (MAX_UNSIGNED_INT ^ ((1 << (L + B)) - 1)))
									 | ((scodon_index & ((1 << L) - 1)) << B)
									 | ((scodon_index >> L) & ((1 << B) - 1));
						scodon[scodon_index] = scodon_buffer;
						scodon_buffer = 0;
					}
					character_index++; 
				}
			}
		}
		fclose(corpus_file_handle);
	}	
	// Calculate statistics.
	sequence_count = sequence_index + 1;
	character_count = character_index;
	sequence_cumulative_length[sequence_count] = character_count;
	sequence_length[sequence_index] = character_count - sequence_cumulative_length[sequence_index];
	scodon_count = (character_count + 16 - 1) >> 4;
	scodon_index = character_index >> 4;
	if (scodon_index < scodon_count)	// There are some corpus nucleotides in the special codon buffer, flush it.
	{
		scodon_index = (scodon_index & (MAX_UNSIGNED_INT ^ ((1 << (L + B)) - 1)))
					 | ((scodon_index & ((1 << L) - 1)) << B)
					 | ((scodon_index >> L) & ((1 << B) - 1));
		scodon[scodon_index] = scodon_buffer; // Now the last special codon might have zeros in its most significant bits. Don't treat such zeros as 'A's.
	}
	scodon_count = ((scodon_count + (1 << (L + B)) - 1) >> (L + B)) << (L + B);	// Now scodon_count is a multiple of 2^(L + B).
}

/**
 * Calculate thread blocks to sequences mapping, i.e. fill block_to_sequence[].
 */
void buildBlockToSequenceMapping()
{
	unsigned int character_index = 0;
	unsigned int sequence_index = 0;
	block_count = scodon_count >> (L + B);	// Calculate the number of thread blocks.	
	for (unsigned int block_index = 0; block_index < block_count; block_index++)
	{
		while (character_index >= sequence_cumulative_length[sequence_index + 1]) sequence_index++;
		block_to_sequence[block_index] = sequence_index;
		character_index += (1 << (L + B + 4)); // One thread block processes 1 << (L + B) special codons, and each special codon encodes 1 << 4 characters.
	}
}

/**
 * Parse the query file. It is the second argument to the program. By default it is Query.txt.
 * Each line of the query file represents a query, which consists of a pattern followed by an edit distance.
 * @param[in] query_file The query file specifying all the queries.
 */
void parseQueryFile(const char *query_file)	// Save the queries into pattern[] and edit_distance[]. Save query_count.
{
	FILE *query_file_handle;	// The handle to query file.
	unsigned int query_length;	// The length of each query.
	query_count = 0;
	query_file_handle = fopen(query_file, "r");
	while (fgets(pattern[query_count], MAX_PATTERN_LENGTH + 3, query_file_handle) != NULL)
	{
		query_length = strlen(pattern[query_count]) - 1; // Exclude the last character which is '\n' of ASCII code '\10'.
		edit_distance[query_count] = pattern[query_count][query_length - 1] - 48;	// char '0' is '\48'.
		pattern[query_count][query_length - 1] = '\0';
		query_count++;
	}
	fclose(query_file_handle);
}

/**
 * The entry of the CUDAagrep program.
 * @param[in] argc Number of arguments to the program. It should be equal to 5. Command: CUDAagrep Corpus.txt Query.txt Result.csv Log.csv
 * @param[in] argv Arguments to the program.
 */
int main(int argc, char** argv)
{
	// About query.
	const char *current_pattern;	// The pattern of current query.
	unsigned int m;	// Pattern length.
	unsigned int k;	// Edit distance.

	// About agrep algorithm.
	unsigned int       mask_array_32[CHARACTER_CARDINALITY];	// The 32-bit mask array of pattern.
	unsigned long long mask_array_64[CHARACTER_CARDINALITY];	// The 64-bit mask array of pattern.
	unsigned int       test_bit_32;	// The test bit for determining matches of patterns of length 32.
	unsigned long long test_bit_64;	// The test bit for determining matches of patterns of length 64.

	// About CUDA implementation.
	unsigned int *scodon_device;	// CUDA global memory pointer pointing to the special codon array.
	unsigned int *match_device;	// CUDA global memory pointer pointing to the match array.
	unsigned int match[MAX_MATCH_COUNT];	// The matches returned by the CUDA agrep kernel.
	unsigned int match_block_value;		// The value of the block of the current matching position. Used as match_block_value = matching_block[i];
	unsigned int m_minus_k;	// m_minus_k = m - k; Used to determine whether a match is across two consecutive sequences.
	unsigned int m_plus_k;	// m_plus_k  = m + k; Used to determine whether a match is across two consecutive sequences.

	// About result.
	unsigned int match_sequence [MAX_MATCH_COUNT];	// The sequence  indexes of matches.
	unsigned int match_character[MAX_MATCH_COUNT];	// The character indexes of matches, i.e. ending positions of matches.
	unsigned int match_sequence_value;	// The sequence  index of a match.
	unsigned int match_character_value;	// The character index of a match.
	unsigned int match_count;	// Actual number of matches in the match array. match_count <= potential_match_count should always holds.
	unsigned int filtered_match_count;	// Actual number of matches after filtering.
	FILE *result_file_handle;	// The handle to result file "Result.csv".
	FILE *log_file_handle;		// The handle to log file "Log.csv".

	// About timing.
	unsigned int timer;	// Each timer has a unique id.
	float query_time;	// The querying time of one single query.

	printf("CUDAagrep v1.3 (Build Januar 13, 2011)\n");
	printf("Hongjian Li (hjli@cse.cuhk.edu.hk)\n");
	if (argc == 1) // argc should be equal to 5 in order to collect all the required arguments.
	{
		printf("Usage: CUDAagrep [InputCorpus] [InputQuery] [OutputResult] [OutputLog]\n");
		printf("e.g.:  CUDAagrep Corpus.txt Query.txt Result.csv Log.csv\n");
		return 0;
	}

	cutilCheckError(cutCreateTimer(&timer));	// Create a timer.
	printf("Reading and parsing corpus files and queries ... ");
	cutilCheckError(cutStartTimer(timer));	// Start timing.
	parseCorpusContainerFile(argv[1]);	// Parse the corpus container file. It is the first argument to the program. By default it is Corpus.txt.
	buildSpecialCodonArray();	// Build the special codon array. Corpus files are enumerated. Statistics of the corpus are calculated.
	buildBlockToSequenceMapping();	// Calculate thread block to sequence index mapping, i.e. fill block_to_sequence[].
	parseQueryFile(argv[2]);	// Parse the query file. It is the second argument to the program. By default it is Query.txt.
	cutilCheckError(cutStopTimer(timer));	// Stop timing.
	printf("done. %u queries. (%.2f ms)\n", query_count, cutGetTimerValue(timer));

	// Allocate the special codon array.
	cutilSafeCall(cudaMalloc((void**)&scodon_device, sizeof(unsigned int) * scodon_count));
	// Transfer the special codon array from host to device.
	cutilSafeCall(cudaMemcpy(scodon_device, scodon, sizeof(unsigned int) * scodon_count, cudaMemcpyHostToDevice));
	// Allocate the match array.
	cutilSafeCall(cudaMalloc((void**)&match_device, sizeof(unsigned int) * MAX_MATCH_COUNT));
	// Transfer necessary arguments to CUDA constant memory.
	initAgrepKernel(scodon_device, character_count, match_device, MAX_MATCH_COUNT);

	// Open result file and log file for writing headers.
	result_file_handle = fopen(argv[3], "w");
	fprintf(result_file_handle, "Query Index,Match Index,Sequence Index,Ending Position,Sequence Header\n");
	log_file_handle = fopen(argv[4], "w");
	fprintf(log_file_handle,    "Query Index,Pattern,Edit Distance,Number of Matches,Querying Time (ms)\n");

	// Use query_index to enumerate pattern[] and edit_distance[].
	for (unsigned int query_index = 0; query_index < query_count; query_index++)
	{
		current_pattern = pattern[query_index];	// Given a pattern, three variables can be derived: m, mask_array, matching_bit.
		m = strlen(current_pattern);	// Pattern length.
		k = edit_distance[query_index];	// Edit distance.
		printf("Query %u: searching for %s with edit distance %u ... ", query_index, current_pattern, k);
		cutilCheckError(cutResetTimer(timer));	// Reset timer.
		cutilCheckError(cutStartTimer(timer));	// Start timing.
		if (m <= 32)
		{
			// Derive mask_array and test_bit.
			memset(mask_array_32, 0, sizeof(unsigned int) * CHARACTER_CARDINALITY);
			for (unsigned int i = 0; i < m; i++)	// Derive the mask array of current pattern.
			{
				unsigned int j = (unsigned int)1 << i;
				if (current_pattern[i] == 'N')
				{
					mask_array_32[0] |= j;
					mask_array_32[1] |= j;
					mask_array_32[2] |= j;
					mask_array_32[3] |= j;
				}
				else
				{
					mask_array_32[encode(current_pattern[i])] |= j;
				}
			}
			mask_array_32[0] ^= MAX_UNSIGNED_INT;
			mask_array_32[1] ^= MAX_UNSIGNED_INT;
			mask_array_32[2] ^= MAX_UNSIGNED_INT;
			mask_array_32[3] ^= MAX_UNSIGNED_INT;
			test_bit_32 = (unsigned int)1 << (m - 1);	// The test bit for determining matches of patterns of length 32.
			transferMaskArray32(mask_array_32, test_bit_32);
		}
		else // m > 32
		{
			// Derive mask_array and test_bit.
			memset(mask_array_64, 0, sizeof(unsigned long long) * CHARACTER_CARDINALITY);
			for (unsigned int i = 0; i < m; i++)	// Derive the mask array of current pattern.
			{
				unsigned long long j = (unsigned long long)1 << i;
				if (current_pattern[i] == 'N')
				{
					mask_array_64[0] |= j;
					mask_array_64[1] |= j;
					mask_array_64[2] |= j;
					mask_array_64[3] |= j;
				}
				else
				{
					mask_array_64[encode(current_pattern[i])] |= j;
				}
			}
			mask_array_64[0] ^= MAX_UNSIGNED_LONG_LONG;
			mask_array_64[1] ^= MAX_UNSIGNED_LONG_LONG;
			mask_array_64[2] ^= MAX_UNSIGNED_LONG_LONG;
			mask_array_64[3] ^= MAX_UNSIGNED_LONG_LONG;
			test_bit_64 = (unsigned long long)1 << (m - 1);	// The test bit for determining matches of patterns of length 64.
			transferMaskArray64(mask_array_64, test_bit_64);
		}
		invokeAgrepKernel(m, k, block_count);
		cutilCheckMsg("CUDA agrep kernel execution failed.");
		// Don't waste the CPU time before waiting for the CUDA agrep kernel to exit.
		m_minus_k = m - k;
		m_plus_k  = m + k;
		filtered_match_count = 0;
		cutilSafeCall(cudaThreadSynchronize());	// Block until the CUDA agrep kernel completes.
		getMatchCount(&match_count);
		if (match_count > MAX_MATCH_COUNT) match_count = MAX_MATCH_COUNT;	// If the number of matches exceeds MAX_MATCH_COUNT, only the first MAX_MATCH_COUNT matches will be saved into the result file.
		cutilSafeCall(cudaMemcpy(match, match_device, sizeof(unsigned int) * match_count, cudaMemcpyDeviceToHost));
		// Enumerate each match.
		for (unsigned int i = 0; i < match_count; i++)
		{
			// match[i] is absolute. Decompose it into relative match_sequence_value and relative match_character_value.
			match_character_value = match[i];	// The absolute ending position of current match.
			match_block_value = match_character_value >> (L + B + 4);	// Derive the thread block that finds the current match.
			match_sequence_value = block_to_sequence[match_block_value];	// Use block-to-sequence mapping to get the nearest sequence index.
			while (match_character_value >= sequence_cumulative_length[match_sequence_value + 1]) match_sequence_value++; // Now match_sequence_value is the sequence index of match[i].
			match_character_value -= sequence_cumulative_length[match_sequence_value];	// Now match_character_value is the character index of match[i].

			if (match_character_value + 1 < m_minus_k) continue; // The current match must be across two consecutive sequences. It is thus an invalid matching.
			match_sequence [filtered_match_count] = match_sequence_value;
			match_character[filtered_match_count] = match_character_value;
			filtered_match_count++;
		}
		cutilCheckError(cutStopTimer(timer));	// Stop timing.
		query_time = cutGetTimerValue(timer);	// Obtain the querying time of current query.
		printf("done. %u matches found. (%.2f ms)\n", filtered_match_count, query_time);
		fprintf(log_file_handle, "%u,%s,%u,%u,%.2f\n", query_index, current_pattern, k, match_count, query_time);
		// Output each filtered match.
		for (unsigned int i = 0; i < filtered_match_count; i++)
		{
			match_sequence_value  = match_sequence [i];
			match_character_value = match_character[i];
			fprintf(result_file_handle, "%u,%u,%u,%u,%s", query_index, i, match_sequence_value, match_character_value, sequence_header[match_sequence_value]);
			if ((match_sequence_value > 0) && (match_character_value + 1 < m_plus_k)) fprintf(result_file_handle, ",?"); // This match may possibly be across two consecutive sequences.
			fprintf(result_file_handle, "\n");
		}
	}
	// All queries are done. Release resources.
	fclose(log_file_handle);
	fclose(result_file_handle);
	cutilSafeCall(cudaFree(scodon_device));
	cutilSafeCall(cudaFree(match_device));
	cutilCheckError(cutDeleteTimer(timer));
	return 0;
}